BWA -- bwa-mem2

1.BWA简介

BWA(Burrows-Wheeler Aligner)主要是将reads比对到大型基因组上，主要功能是：序列比对。
首先通过BWT(Burrows-Wheeler Transformation，BWT压缩算法)为大型参考基因组建立索引，然后将reads比对到基因组。特点是快速、准确、省内存。
由三种类似算法组成：BWA-backtrack，BWA-SW和BWA-MEM。
首推BWA-MEM。

1.1三种算法的适用范围

1.2BWA参数

序列比对BWA之参数：index, mem, aln, samse, sampe, bwasw

-   bwa index ref.fa # 首先建立基因组索引
-   bwa mem ref.fa reads.fq > aln-se.sam # 调用BWA-MEM
-   bwa mem ref.fa read1.fq read2.fq > aln-pe.sam # 调用BWA-MEM
-   bwa aln ref.fa short_read.fq > aln_sa.sai # 调用BWA-backtrack
-   bwa samse ref.fa aln_sa.sai short_read.fq > aln-se.sam # 调用BWA-backtrack
-   bwa sampe ref.fa aln_sa1.sai aln_sa2.sai read1.fq read2.fq > aln-pe.sam # 调用BWA-backtrack
-   bwa bwasw ref.fa long_read.fq > aln.sam # 调用BWA-SW

-   注意：BWA输入的是fastq/fq的原始测序数据。

1.3Published Articles:

2.bwa-mem2

2.1安装

$ curl -L https://github.com/bwa-mem2/bwa-mem2/releases/download/v2.2.1/bwa-mem2-2.2.1_x64-linux.tar.bz2 | tar jxf -
$ bwa-mem2-2.2.1_x64-linux/bwa-mem2 index ref.fa
$ bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem ref.fa read1.fq read2.fq > out.sam

2.2使用

# 建立索引
$ bwa-mem2-2.2.1_x64-linux/bwa-mem2 index ref.fa

# bwa比对
$ nohup time bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -t 24 ref.fa JLTG01.filter.R1.fq.gz JLTG01.filter.R2.fq.gz 1>JLTG01.bwa.sam 2>JLTG01.bwa.log &

-   -R "@RG\tID:W2018001\tPL:ILLUMINA\tLB:W2018001\tSM:W2018001" # 文件头信息，sam文件的，PL是有要求的，其余的自己设就好了

3.批量运行与分步详解

以下内容整理自实际项目经验，补充批量操作和详细注释。

3.1批量运行实例

以GP-20221121-5291项目5个猴子样品为例

Pasted image 20230510172857.png

$ for i in {J5,Z3,Z4,ZH1,ZH2};do echo nohup time /share/nas1/yuj/software/bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -t 24 GCF_003339765.1_Mmul_10_genomic.fna  /share/nas1/seqloader/xianliti/GP-20221121-5291_zhongguoyixuekexueyuan_35samples_hou_xianliti/2.clean_data/$i/*_R1_000.fastq.gz /share/nas1/seqloader/xianliti/GP-20221121-5291_zhongguoyixuekexueyuan_35samples_hou_xianliti/2.clean_data/$i/*_R2_000.fastq.gz  "1>"$i/$i.bwa.sam  "2>"$i.bwa.log "&&" samtools flagstat -@ 10 $i/$i.bwa.sam ">" $i.bwa.flagstat.log "&";done
# 这条命令效果是显示上图蓝色的5条命令,一条命令对应一个样本,可以一次性复制回车运行
# {J5,Z3,Z4,ZH1,ZH2}是样品名
# 具体解释见下面的分步讲解

GP-20230306-5682项目11个样挑了5个样

$ for i in {DCYD01,JLWXH01,KDWN01,LP1,DQBST201};do echo $i;nohup time bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -t 10 Primula_veris_subsp._veris.faa /share/nas1/seqloader/ddrad/GP-20230306-5682_huaibeishifandaxue_112samples_baochunhua_ddrad/all_data/2.clean_data/$i.filter.R1.fq.gz /share/nas1/seqloader/ddrad/GP-20230306-5682_huaibeishifandaxue_112samples_baochunhua_ddrad/all_data/2.clean_data/$i.filter.R2.fq.gz 1>$i.bwa.sam 2>$i.bwa.log && samtools flagstat -@ 10 $i.bwa.sam >$i.bwa.flagstat.log & done
# 这条命令的效果是直接运行了生成的5条命令

3.2分步拆解（带注释）

## .gz文件 解压缩
for i in *.gz;do echo $i;gunzip -c $i > ${i%.gz};done

建立索引

$ /share/nas1/yuj/software/bwa-mem2-2.2.1_x64-linux/bwa-mem2  index  Primula_veris_subsp._veris.faa(参考基因组)

bwa比对

$ nohup time /share/nas1/yuj/software/bwa-mem2-2.2.1_x64-linux/bwa-mem2 mem -t 24 Primula_veris_subsp._veris.faa(参考基因组)  JLTG01.filter.R1.fq.gz  JLTG01.filter.R2.fq.gz  1>JLTG01.bwa.sam  2>JLTG01.bwa.log  &

- time这个是为了统计运行时间,可加可不加
- /share/nas1/yuj/software/bwa-mem2-2.2.1_x64-linux/bwa-mem2程序路径
- -t表示线程数
- Primula_veris_subsp._veris.faa是参考基因组
- JLTG01.filter.R1.fq.gz测序reads
- JLTG01.filter.R2.fq.gz另一端测序reads
- 生成JLTG01.bwa.sam是比对上的reads

samtools统计

$ samtools flagstat -@ 10 JLTG01.bwa.sam > JLTG01.bwa.flagstat.log &
- -@ 10是线程数
- JLTG01.bwa.sam上一个程序结果
- JLTG01.bwa.flagstat.log统计后的结果,看第五行

3.3实用小命令

# 显示样品名
$ for i in `ls /share/nas1/seqloader/ddrad/GP-20230306-5682_huaibeishifandaxue_112samples_baochunhua_ddrad/all_data/2.clean_data`;do echo $i;done |awk -F "."  'BEGIN {max = 0} {if ($1 != max) print $1}{max=$1} END {print "done"}'

4.RSeQC

RSeQC是一个功能强大的软件，里面有很多实用的小工具，其中的bam_stat就是一个实用的bam/sam结果统计工具，安装方式也是相当简单了，就是一个python的包，支持python2.x和python3.x，这里我选用python3的pip来安装，因为本人习惯使用python3。

4.1安装

$ pip3 install RSeQC

4.2统计

$ bam_stat.py -i JLTG01.bwa.sam

相关笔记